1 Introducción y definición de objetivos

El análisis de datos en el mercado inmobiliario es algo muy común y que lleva muchísimos años desarrollándose. Con el objetivo de predecir que aspectos influyen principalmente en el precio de las casas en Madrid, hemos seleccionado un dataset (kaggle - Madrid House Price) con viviendas en venta de la capital española.

library(readr)
library(dplyr)
library(tidyr)
library(kableExtra)
library(ggplot2)
library(GGally)
library(gridExtra)
library(cowplot)
library(ggcorrplot)
library(gmodels)
library(caret)
library(ggfortify)
library(scales)
library(class)
library(distances)
library(visreg)
library(rpart)
library(rpart.plot)
library(rattle)
library(randomForest)
library(e1071)
library(kernlab)
library(pROC)
library(cluster)
library(tidyverse)
library(stringr)
library(Metrics)
library(factoextra)
library(NbClust)
paste(R.Version()$version.string)
## [1] "R version 4.1.2 (2021-11-01)"

2 Análisis exploratorio inicial

2.1 Lectura y preparación de los datos

En esta primera etapa cargamos el dataset descargado y hacemos una observación superficial de las variables recogidas y sus características. Se representan también las 10 filas iniciales.

mhp <- read_csv("house_price_madrid_14_08_2022.csv")
head(mhp, 10) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
price house_type house_type_2 rooms m2 elevator garage neighborhood district
495000 planta 1 exterior 3 118 TRUE TRUE Chopera Arganzuela
485000 planta 2 exterior 2 82 TRUE TRUE Palos de Moguer Arganzuela
315000 planta 2 exterior 2 72 FALSE FALSE Legazpi Arganzuela
585000 planta 4 exterior 2 174 TRUE TRUE Palos de Moguer Arganzuela
255000 bajo exterior 3 75 FALSE FALSE Acacias Arganzuela
299000 planta 1 exterior 1 69 TRUE FALSE Chopera Arganzuela
265000 planta 3 exterior 2 54 TRUE FALSE Delicias Arganzuela
290000 planta 3 interior 4 69 TRUE FALSE Chopera Arganzuela
660220 planta 2 exterior 3 129 TRUE TRUE Imperial Arganzuela
525000 planta 4 exterior 4 111 TRUE FALSE Chopera Arganzuela
summary(mhp)
##      price           house_type        house_type_2           rooms       
##  Min.   :     725   Length:15975       Length:15975       Min.   : 1.000  
##  1st Qu.:  195000   Class :character   Class :character   1st Qu.: 2.000  
##  Median :  359973   Mode  :character   Mode  :character   Median : 3.000  
##  Mean   :  624233                                         Mean   : 2.848  
##  3rd Qu.:  749000                                         3rd Qu.: 3.000  
##  Max.   :13950000                                         Max.   :41.000  
##        m2         elevator         garage        neighborhood      
##  Min.   :  1.0   Mode :logical   Mode :logical   Length:15975      
##  1st Qu.: 66.0   FALSE:4597      FALSE:11528     Class :character  
##  Median : 93.0   TRUE :11378     TRUE :4447      Mode  :character  
##  Mean   :124.8                                                     
##  3rd Qu.:142.0                                                     
##  Max.   :989.0                                                     
##    district        
##  Length:15975      
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
str(mhp) 
## spc_tbl_ [15,975 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ price       : num [1:15975] 495000 485000 315000 585000 255000 ...
##  $ house_type  : chr [1:15975] "planta 1" "planta 2" "planta 2" "planta 4" ...
##  $ house_type_2: chr [1:15975] "exterior" "exterior" "exterior" "exterior" ...
##  $ rooms       : num [1:15975] 3 2 2 2 3 1 2 4 3 4 ...
##  $ m2          : num [1:15975] 118 82 72 174 75 69 54 69 129 111 ...
##  $ elevator    : logi [1:15975] TRUE TRUE FALSE TRUE FALSE TRUE ...
##  $ garage      : logi [1:15975] TRUE TRUE FALSE TRUE FALSE FALSE ...
##  $ neighborhood: chr [1:15975] "Chopera" "Palos de Moguer" "Legazpi" "Palos de Moguer" ...
##  $ district    : chr [1:15975] "Arganzuela" "Arganzuela" "Arganzuela" "Arganzuela" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   price = col_double(),
##   ..   house_type = col_character(),
##   ..   house_type_2 = col_character(),
##   ..   rooms = col_double(),
##   ..   m2 = col_double(),
##   ..   elevator = col_logical(),
##   ..   garage = col_logical(),
##   ..   neighborhood = col_character(),
##   ..   district = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>

El dataset contiene 15.975 observaciones (correspondientes a una vivienda cada una) y 9 variables (de las cuales son 6 cualitativas y 3 cuantitativas).

A continuación, la descripción de cada una de las variables:

  • price: precio

  • house_type: tipo de vivienda (casa, chalet, piso…)

  • house_type_2: si es exterior o interior

  • rooms: número de habitaciones

  • m2: metros cuadrados

  • elevator: si tiene ascensor

  • garage: si incluye garaje

  • neighborhood: barrio de Madrid

  • district: distrito de Madrid

Tras esta primera observación se harán algunas modificaciones en el dataset para transformar a variables de tipo categórico las que correspondan, así como transformar a binarias (1/0) aquellas con sólo dos valores posibles.

# Preparación de los datos. Convertimos a factores las variables que son categóricas y a 1-0 las binarias.

mhp <- mhp %>% rename(exterior = house_type_2)
mhp$exterior = ifelse(mhp$exterior == "exterior", 1, 0)
mhp$exterior = factor(mhp$exterior, levels = c(1,0))

mhp$elevator = ifelse(mhp$elevator == "TRUE", 1, 0)
mhp$elevator = factor(mhp$elevator, levels = c(1,0))

mhp$garage = ifelse(mhp$garage == "TRUE", 1, 0)
mhp$garage = factor(mhp$garage, levels = c(1,0))

mhp$house_type <- as.factor(mhp$house_type)
mhp$exterior <- as.factor(mhp$exterior)
mhp$elevator <- as.factor(mhp$elevator)
mhp$garage <- as.factor(mhp$garage)
mhp$neighborhood <- as.factor(mhp$neighborhood)
mhp$district <- as.factor(mhp$district)

2.2 Tratamiento de datos faltantes

Convertimos en NA aquellas observaciones que claramente contienen algún error:

  • La casa de 41 habitaciones (probablemente fueran 4 y esté mal imputado)

  • Las que tienen menos de 3 m2 (probablemente mal imputados al usar el punto para separar los miles)

  • La casa con un precio de 725 (claramente equivocado).

mhp$rooms[mhp$rooms > 40] <- NA
mhp$m2[mhp$m2 < 3] <- NA
mhp$price[mhp$price < 1000] <- NA

Por otro lado, la variable house_type_2 tiene 469 filas sin datos. Junto a los NA añadidos en el paso anterior son el 3.2% del total, por lo que optaremos por eliminar estas filas del dataset.

sum(is.na(mhp))/nrow(mhp)
## [1] 0.03205008
mhp <- na.omit(mhp)

2.3 División del dataset

Para finalizar la preparación de los datos, dividiremos nuestro dataset en tres partes, una de entrenamiento (train), con el que trabajaremos durante todo el análisis, otra de prueba (test) que nos servirá para comprobar la eficacia de los diferentes modelos predictivos que se apliquen, y finalmente uno de validación (validation), que se mantendrá sin usar hasta el último día del proyecto y con el cual validaremos el modelo final.

set.seed(108)
numero_total = nrow(mhp)
# Porcentajes de train, test y validation
w_train = .5
w_test = .25
w_validation = 1 - (w_train + w_test)

# Todos los índices
indices = seq(1:numero_total) 

# Muestreo
indices_train = sample(1:numero_total, numero_total * w_train)
indices_test = sample(indices[-indices_train], numero_total * w_test)
indices_validation = indices[-c(indices_train,indices_test)]

# Agrupamos

mhp_train = mhp[indices_train,]
mhp_test = mhp[indices_test,]
mhp_validation = mhp[indices_validation,]

3 Análisis univariante

3.1 Análisis de variables cualitativas

Una vez preparado el dataset comenzamos con el análisis de las variables cualitativas, que son las siguientes:

Tipo de vivienda

merge(setNames(as.data.frame(table(mhp_train$house_type)), c("house_type", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$house_type))*100, 2)), c("house_type", "prop (%)"))
) %>%
  arrange(desc(count)) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
house_type count prop (%)
planta 1 1581 20.45
planta 2 1332 17.23
planta 3 1180 15.26
bajo 996 12.88
planta 4 801 10.36
planta 5 545 7.05
planta 6 317 4.10
planta 7 196 2.53
chalet 180 2.33
entreplanta 131 1.69
planta 8 119 1.54
planta 9 78 1.01
casa 65 0.84
semi-sotano 48 0.62
planta 10 44 0.57
planta 12 23 0.30
planta 13 22 0.28
planta 11 21 0.27
sotano 19 0.25
planta 14 10 0.13
planta -1 8 0.10
planta 15 6 0.08
planta 16 4 0.05
planta 18 2 0.03
planta 19 2 0.03
planta 20 2 0.03
ggplot(mhp_train, aes(house_type)) +
  geom_bar(fill = "#0BB363") +
  coord_flip() +
  labs(x = "Tipo de vivienda", y = "Número de viviendas", title = "Viviendas por tipo") +
  theme(plot.title = element_text(hjust = 0.5))

Como era de esperar, las viviendas más habituales son pisos en plantas no demasiado altas (bajos, primeros, segundos…). El número de casas y chalets también es pequeño en proporción.

Orientación exterior

merge(setNames(as.data.frame(table(mhp_train$exterior)), c("exterior", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$exterior))*100, 2)), c("exterior", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
exterior count prop (%)
0 848 10.97
1 6884 89.03
ggplot(mhp_train, aes(exterior)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Exterior", y = "Número de vivienda de viviendas", title = "Viviendas exteriores") +
  theme(plot.title = element_text(hjust = 0.5))

Son una inmensa mayoría los pisos que dan a exterior.

Ascensor

merge(setNames(as.data.frame(table(mhp_train$elevator)), c("elevator", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$elevator))*100, 2)), c("elevator", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
elevator count prop (%)
0 2169 28.05
1 5563 71.95
ggplot(mhp_train, aes(elevator)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Ascensor", y = "Número de viviendas", title = "Viviendas con ascensor") +
  theme(plot.title = element_text(hjust = 0.5))

También es más frecuente que tengan ascensor, aunque casi el 30% no lo tienen.

Garaje

merge(setNames(as.data.frame(table(mhp_train$garage)), c("garage", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$garage))*100, 2)), c("garage", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
garage count prop (%)
0 5570 72.04
1 2162 27.96
ggplot(mhp_train, aes(garage)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Garaje", y = "Número de viviendas", title = "Viviendas con garaje") +
  theme(plot.title = element_text(hjust = 0.5))

Parece normal que la mayoría de pisos no tuvieran garaje incluido.

Barrio

Esta variable no será muy útil para el análisis, ya que en muchos casos son descripciones de la vivienda hechas por el usuario, por lo que habría que hacer un tratamiento previo para poder agruparlas. Algo innecesario ya que ese mismo efecto se consigue con el siguiente campo: distritos.

merge(setNames(as.data.frame(table(mhp_train$neighborhood)), c("neighborhood", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$neighborhood))*100, 2)), c("neighborhood", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
neighborhood count prop (%)
12 de Octubre-Orcasur 19 0.25
Abrantes 54 0.70
Acacias 50 0.65
Adelfas 30 0.39
Aeropuerto 2 0.03
Águilas 33 0.43
Alameda de Osuna 32 0.41
Almagro 140 1.81
Almendrales 42 0.54
Aluche 56 0.72
Ambroz 30 0.39
Amposta 18 0.23
Apóstol Santiago 19 0.25
Arapiles 67 0.87
Aravaca 47 0.61
Arcos 26 0.34
Argüelles 142 1.84
Arroyo del Fresno 11 0.14
Atalaya 6 0.08
Ático en Almagro 8 0.10
Ático en Arapiles 1 0.01
Ático en Aravaca 0 0.00
Ático en Arcos 1 0.01
Ático en Argüelles 0 0.00
Ático en Bernabéu-Hispanoamérica 1 0.01
Ático en Campo de las Naciones-Corralejos 1 0.01
Ático en Castellana 10 0.13
Ático en Castilla 4 0.05
Ático en Chueca-Justicia 1 0.01
Ático en Ciudad Universitaria 1 0.01
Ático en Colina 1 0.01
Ático en Concepción 1 0.01
Ático en Conde Orgaz-Piovera 1 0.01
Ático en Costillares 1 0.01
Ático en Cuatro Caminos 0 0.00
Ático en Delicias 1 0.01
Ático en El Cañaveral 1 0.01
Ático en El Viso 8 0.10
Ático en Fuente del Berro 0 0.00
Ático en Gaztambide 2 0.03
Ático en Goya 3 0.04
Ático en Guindalera 2 0.03
Ático en Horcajo 1 0.01
Ático en Ibiza 2 0.03
Ático en Jerónimos 1 0.01
Ático en Las Tablas 1 0.01
Ático en Lista 1 0.01
Ático en Malasaña-Universidad 1 0.01
Ático en Niño Jesús 1 0.01
Ático en Nueva España 4 0.05
Ático en Nuevos Ministerios-Ríos Rosas 3 0.04
Ático en Pacífico 1 0.01
Ático en Palacio 0 0.00
Ático en Palomas 1 0.01
Ático en Peñagrande 1 0.01
Ático en Pinar del Rey 0 0.00
Ático en Pueblo Nuevo 1 0.01
Ático en Puerta Bonita 0 0.00
Ático en Puerta del Ángel 1 0.01
Ático en Recoletos 3 0.04
Ático en Rejas 0 0.00
Ático en Rosas 1 0.01
Ático en Salvador 1 0.01
Ático en San Juan Bautista 1 0.01
Ático en San Pascual 1 0.01
Ático en Simancas 0 0.00
Ático en Sol 3 0.04
Ático en Trafalgar 1 0.01
Ático en Valdeacederas 1 0.01
Ático en Valdebebas - Valdefuentes 1 0.01
Ático en Valdebernardo - Valderrivas 0 0.00
Ático en Valdemarín 2 0.03
Ático en Vallehermoso 2 0.03
Bellas Vistas 89 1.15
Bernabéu-Hispanoamérica 91 1.18
Berruguete 82 1.06
Buena Vista 40 0.52
Butarque 14 0.18
Campamento 16 0.21
Campo de las Naciones-Corralejos 16 0.21
Canillas 43 0.56
Canillejas 78 1.01
Casa de Campo 18 0.23
Casa o chalet independiente en Aravaca 7 0.09
Casa o chalet independiente en Bernabéu-Hispanoamérica 2 0.03
Casa o chalet independiente en Canillas 2 0.03
Casa o chalet independiente en Canillejas 1 0.01
Casa o chalet independiente en Casa de Campo 0 0.00
Casa o chalet independiente en Casco Histórico de Barajas 0 0.00
Casa o chalet independiente en Castilla 0 0.00
Casa o chalet independiente en Ciudad Universitaria 5 0.06
Casa o chalet independiente en Conde Orgaz-Piovera 7 0.09
Casa o chalet independiente en El Plantío 1 0.01
Casa o chalet independiente en El Viso 4 0.05
Casa o chalet independiente en Fuentelarreina 0 0.00
Casa o chalet independiente en Guindalera 0 0.00
Casa o chalet independiente en Mirasierra 4 0.05
Casa o chalet independiente en Niño Jesús 1 0.01
Casa o chalet independiente en Nueva España 6 0.08
Casa o chalet independiente en Opañel 0 0.00
Casa o chalet independiente en Orcasitas 0 0.00
Casa o chalet independiente en Peñagrande 4 0.05
Casa o chalet independiente en Rejas 1 0.01
Casa o chalet independiente en Tres Olivos - Valverde 2 0.03
Casa o chalet independiente en Valdemarín 7 0.09
Casa o chalet independiente en Valdezarza 0 0.00
Casa o chalet independiente en Villaverde Alto 0 0.00
Casa o chalet independiente en Vista Alegre 0 0.00
Casco Histórico de Barajas 22 0.28
Casco Histórico de Vallecas 87 1.13
Casco Histórico de Vicálvaro 35 0.45
Castellana 172 2.22
Castilla 62 0.80
Chalet adosado en Alameda de Osuna 1 0.01
Chalet adosado en Aravaca 0 0.00
Chalet adosado en Arroyo del Fresno 1 0.01
Chalet adosado en Bernabéu-Hispanoamérica 1 0.01
Chalet adosado en Buena Vista 1 0.01
Chalet adosado en Campo de las Naciones-Corralejos 1 0.01
Chalet adosado en Canillas 1 0.01
Chalet adosado en Canillejas 0 0.00
Chalet adosado en Ciudad Jardín 1 0.01
Chalet adosado en Ciudad Universitaria 0 0.00
Chalet adosado en Concepción 1 0.01
Chalet adosado en Conde Orgaz-Piovera 5 0.06
Chalet adosado en El Plantío 2 0.03
Chalet adosado en El Viso 4 0.05
Chalet adosado en Fuente del Berro 1 0.01
Chalet adosado en Guindalera 0 0.00
Chalet adosado en Los Cármenes 0 0.00
Chalet adosado en Mirasierra 1 0.01
Chalet adosado en Moscardó 1 0.01
Chalet adosado en Nueva España 3 0.04
Chalet adosado en Numancia 0 0.00
Chalet adosado en Orcasitas 2 0.03
Chalet adosado en Palomas 1 0.01
Chalet adosado en Peñagrande 1 0.01
Chalet adosado en Pinar del Rey 1 0.01
Chalet adosado en Puerta Bonita 0 0.00
Chalet adosado en Rejas 1 0.01
Chalet adosado en Salvador 0 0.00
Chalet adosado en Tres Olivos - Valverde 0 0.00
Chalet adosado en Valdeacederas 0 0.00
Chalet adosado en Valdebebas - Valdefuentes 1 0.01
Chalet adosado en Valdezarza 1 0.01
Chalet en Águilas 0 0.00
Chalet en Aravaca 1 0.01
Chalet en Ciudad Universitaria 1 0.01
Chalet en El Plantío 1 0.01
Chalet en El Viso 0 0.00
Chalet en Entrevías 1 0.01
Chalet en Fuentelarreina 0 0.00
Chalet en Mirasierra 1 0.01
Chalet en Rejas 0 0.00
Chalet en Valdemarín 0 0.00
Chalet pareado en Abrantes 1 0.01
Chalet pareado en Alameda de Osuna 0 0.00
Chalet pareado en Apóstol Santiago 1 0.01
Chalet pareado en Aravaca 3 0.04
Chalet pareado en Arroyo del Fresno 1 0.01
Chalet pareado en Bernabéu-Hispanoamérica 0 0.00
Chalet pareado en Campo de las Naciones-Corralejos 0 0.00
Chalet pareado en Canillas 0 0.00
Chalet pareado en Canillejas 0 0.00
Chalet pareado en Ciudad Jardín 0 0.00
Chalet pareado en Ciudad Universitaria 1 0.01
Chalet pareado en Colina 1 0.01
Chalet pareado en Conde Orgaz-Piovera 4 0.05
Chalet pareado en El Plantío 2 0.03
Chalet pareado en El Viso 3 0.04
Chalet pareado en Fuentelarreina 0 0.00
Chalet pareado en Guindalera 0 0.00
Chalet pareado en Mirasierra 3 0.04
Chalet pareado en Nueva España 0 0.00
Chalet pareado en Palomas 4 0.05
Chalet pareado en Peñagrande 1 0.01
Chalet pareado en Rejas 1 0.01
Chalet pareado en Tres Olivos - Valverde 0 0.00
Chalet pareado en Valdebebas - Valdefuentes 1 0.01
Chopera 34 0.44
Chueca-Justicia 67 0.87
Ciudad Jardín 47 0.61
Ciudad Universitaria 40 0.52
Colina 25 0.32
Comillas 52 0.67
Concepción 49 0.63
Conde Orgaz-Piovera 35 0.45
Costillares 53 0.69
Cuatro Caminos 98 1.27
Cuatro Vientos 1 0.01
Cuzco-Castillejos 98 1.27
Delicias 58 0.75
Dúplex en Almagro 6 0.08
Dúplex en Apóstol Santiago 0 0.00
Dúplex en Arapiles 0 0.00
Dúplex en Aravaca 3 0.04
Dúplex en Arcos 0 0.00
Dúplex en Argüelles 0 0.00
Dúplex en Arroyo del Fresno 1 0.01
Dúplex en Bernabéu-Hispanoamérica 0 0.00
Dúplex en Berruguete 1 0.01
Dúplex en Butarque 0 0.00
Dúplex en Canillas 1 0.01
Dúplex en Canillejas 1 0.01
Dúplex en Casco Histórico de Barajas 0 0.00
Dúplex en Casco Histórico de Vallecas 0 0.00
Dúplex en Castellana 1 0.01
Dúplex en Castilla 1 0.01
Dúplex en Chueca-Justicia 1 0.01
Dúplex en Colina 1 0.01
Dúplex en Concepción 2 0.03
Dúplex en Conde Orgaz-Piovera 1 0.01
Dúplex en Costillares 0 0.00
Dúplex en Cuzco-Castillejos 1 0.01
Dúplex en El Viso 4 0.05
Dúplex en Entrevías 1 0.01
Dúplex en Fontarrón 1 0.01
Dúplex en Fuente del Berro 0 0.00
Dúplex en Goya 1 0.01
Dúplex en Imperial 1 0.01
Dúplex en Jerónimos 4 0.05
Dúplex en Las Tablas 2 0.03
Dúplex en Legazpi 2 0.03
Dúplex en Lista 0 0.00
Dúplex en Malasaña-Universidad 0 0.00
Dúplex en Mirasierra 1 0.01
Dúplex en Moscardó 1 0.01
Dúplex en Nueva España 8 0.10
Dúplex en Nuevos Ministerios-Ríos Rosas 1 0.01
Dúplex en Numancia 2 0.03
Dúplex en Opañel 1 0.01
Dúplex en Pacífico 3 0.04
Dúplex en Palacio 1 0.01
Dúplex en Pinar del Rey 0 0.00
Dúplex en Prosperidad 1 0.01
Dúplex en Puerta Bonita 0 0.00
Dúplex en Puerta del Ángel 0 0.00
Dúplex en Quintana 0 0.00
Dúplex en Recoletos 3 0.04
Dúplex en Rejas 0 0.00
Dúplex en San Isidro 1 0.01
Dúplex en San Juan Bautista 6 0.08
Dúplex en Sanchinarro 1 0.01
Dúplex en Simancas 2 0.03
Dúplex en Trafalgar 2 0.03
Dúplex en Tres Olivos - Valverde 1 0.01
Dúplex en Valdeacederas 3 0.04
Dúplex en Valdebernardo - Valderrivas 0 0.00
Dúplex en Valdemarín 4 0.05
Dúplex en Vallehermoso 1 0.01
Dúplex en Ventilla-Almenara 0 0.00
Dúplex en Villaverde Alto 1 0.01
El Cañaveral 50 0.65
El Pardo 3 0.04
El Plantío 3 0.04
El Viso 109 1.41
Ensanche de Vallecas - La Gavia 60 0.78
Entrevías 38 0.49
Estrella 27 0.35
Fontarrón 31 0.40
Fuente del Berro 45 0.58
Fuentelarreina 9 0.12
Gaztambide 77 1.00
Goya 257 3.32
Guindalera 97 1.25
Hellín 17 0.22
Horcajo 3 0.04
Huertas-Cortes 48 0.62
Ibiza 56 0.72
Imperial 48 0.62
Jerónimos 62 0.80
La Paz 48 0.62
Las Tablas 39 0.50
Lavapiés-Embajadores 115 1.49
Legazpi 38 0.49
Lista 134 1.73
Los Ángeles 43 0.56
Los Berrocales 3 0.04
Los Cármenes 50 0.65
Los Rosales 34 0.44
Lucero 43 0.56
Malasaña-Universidad 126 1.63
Marroquina 17 0.22
Media Legua 20 0.26
Mirasierra 15 0.19
Montecarmelo 14 0.18
Moscardó 47 0.61
Niño Jesús 24 0.31
Nueva España 76 0.98
Nuevos Ministerios-Ríos Rosas 91 1.18
Numancia 90 1.16
Opañel 54 0.70
Orcasitas 38 0.49
Pacífico 83 1.07
Palacio 88 1.14
Palomas 20 0.26
Palomeras Bajas 42 0.54
Palomeras sureste 37 0.48
Palos de Moguer 73 0.94
Pau de Carabanchel 26 0.34
Pavones 6 0.08
Peñagrande 77 1.00
Pilar 63 0.81
Pinar del Rey 67 0.87
Portazgo 43 0.56
Pradolongo 32 0.41
Prosperidad 86 1.11
Pueblo Nuevo 132 1.71
Puerta Bonita 89 1.15
Puerta del Ángel 137 1.77
Quintana 68 0.88
Recoletos 181 2.34
Rejas 59 0.76
Rosas 30 0.39
Salvador 21 0.27
San Cristóbal 30 0.39
San Diego 120 1.55
San Fermín 33 0.43
San Isidro 91 1.18
San Juan Bautista 39 0.50
San Pascual 48 0.62
Sanchinarro 61 0.79
Santa Eugenia 16 0.21
Simancas 60 0.78
Sol 66 0.85
Timón 14 0.18
Trafalgar 97 1.25
Tres Olivos - Valverde 51 0.66
Valdeacederas 104 1.35
Valdebebas - Valdefuentes 42 0.54
Valdebernardo - Valderrivas 22 0.28
Valdemarín 19 0.25
Valdezarza 42 0.54
Vallehermoso 72 0.93
Ventas 103 1.33
Ventilla-Almenara 46 0.59
Villaverde Alto 87 1.13
Vinateros 20 0.26
Virgen del Cortijo - Manoteras 31 0.40
Vista Alegre 102 1.32
Zofío 25 0.32

Distrito

En este caso los datos si están correctamente recogidos, repartiendo las casas entre los 21 distritos de Madrid.Casi todos con una representación considerable, en lo que destaca el barrio de Salamanca, con casi el doble de viviendas que el segundo que más tiene, Chamberí.

merge(setNames(as.data.frame(table(mhp_train$district)), c("district", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$district))*100, 2)), c("district", "prop (%)"))
) %>%
  arrange(desc(count)) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
district count prop (%)
barrio de salamanca 911 11.78
chamberi 571 7.38
ciudad lineal 540 6.98
chamartin 526 6.80
tetuan 523 6.76
centro 517 6.69
carabanchel 512 6.62
puente-de-vallecas 374 4.84
fuencarral 356 4.60
moncloa 353 4.57
hortaleza 352 4.55
latina 337 4.36
san-blas 319 4.13
Arganzuela 305 3.94
retiro 295 3.82
usera 240 3.10
villaverde 209 2.70
villa-de-vallecas 163 2.11
vicalvaro 141 1.82
moratalaz 99 1.28
barajas 89 1.15
ggplot(mhp_train, aes(district)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Distrito", y = "Número de viviendas", title = "Viviendas por distrito") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

3.2 Análisis de variables cuantitativas

Precio

data.frame(summarise(mhp_train,
                     min = min(price),
                     max = max(price),
                     median = median(price),
                     mean = mean(price),
                     sd = sd(price))) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
min max median mean sd
47500 8950000 369000 631920.1 753136.6
ggplot(mhp_train, aes(y = price)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Precio", title = "Boxplot de precios") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(price)) +
  geom_histogram(aes(y = ..count..), bins = 50, position = "dodge", fill = "#0BB363") +
  labs(x = "Precio", y = "Número de viviendas", title = "Viviendas por Precio") +
  scale_x_continuous(breaks = pretty(mhp_train$price, n = 10), labels = comma_format()) +
  scale_y_continuous(labels = function(x) sprintf("%.0f", x)) +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30))

Habitaciones

data.frame(summarise(mhp_train,
                     min = min(rooms),
                     max = max(rooms),
                     median = median(rooms),
                     mean = mean(rooms),
                     sd = sd(rooms))) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
min max median mean sd
1 17 3 2.848552 1.302995
ggplot(mhp_train, aes(y = rooms)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Habitaciones", title = "Boxplot de habitaciones") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(rooms)) +
  geom_histogram(aes(), bins = 16, position = "dodge", fill = "#0BB363") +
  geom_density(alpha=.2, fill = "red") +
  labs(x = "Habitaciones", y = "Número de viviendas", title = "Viviendas por número de habitaciones") +
  theme(plot.title = element_text(hjust = 0.5))

data.frame(summarise(mhp_train,
                     min = min(m2),
                     max = max(m2),
                     median = median(m2),
                     mean = mean(m2),
                     sd = sd(m2))) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
min max median mean sd
20 986 95 127.1208 102.497
ggplot(mhp_train, aes(y = m2)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Metros cuadrados", title = "Boxplot de superficies") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(m2)) +
  geom_histogram(aes(y=..count..), bins = 50, position = "dodge", fill = "#0BB363") +
  geom_density(alpha=.3, fill = "red") +
  labs(x = "Metros cuadrados", y = "Número de viviendas", title = "Viviendas por metros cuadrados") +
  theme(plot.title = element_text(hjust = 0.5))

4 Análisis multivariante

Una vez analizadas todas las variables de forma individual, podemos buscar algunas relaciones entre ellas.

Como se ve en los gráficos inferiores hay una relción fuerte entre el precio, los metros cuadradados y el número de habitaciones.

plot_grid(
  ggcorrplot(cor(mhp_train %>% select_if(is.numeric)),type = "lower", lab=TRUE),
  
  ggplot(mhp_train, aes(x = m2, y = price)) +
  geom_point() +
  geom_smooth() +
  ggtitle('Reación precio y m²') +
  theme(plot.title = element_text(hjust = 0.5)),
  
  ggplot(mhp_train, aes(x = rooms, y = price)) +
  geom_point() +
  geom_smooth() +
  ggtitle('Reación precio y habitaciones') +
  theme(plot.title = element_text(hjust = 0.5)), 
  
  ggplot(mhp_train, aes(x = rooms, y = m2)) +
  geom_point() +
  geom_smooth() +
  ggtitle('Reación m² y habitaciones') +
  theme(plot.title = element_text(hjust = 0.5)), 
  
  nrow = 2
)

Introduciendo algunas de las variables categóricas se puede observar como varía el precio con la superficie o el número de habitaciones dependiendo de si incluye garaje, ascensor u orientación exterior.

plot_grid(
  ggplot(mhp_train, aes(x = m2, y = price, colour = exterior)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-superficie, exterior') ,
  
  ggplot(mhp_train, aes(x = rooms, y = price, colour = exterior)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) + 
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-habitaciones, exterior') ,
  
  nrow =1
)

plot_grid(
  ggplot(mhp_train, aes(x = m2, y = price, colour = elevator)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-superficie, ascensor'),
  
  ggplot(mhp_train, aes(x = rooms, y = price, colour = elevator)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-habitación, ascensor'),
  
  nrow =1
)

plot_grid(
  ggplot(mhp_train, aes(x = m2, y = price, colour = garage)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-superficie, garaje'),
  
  ggplot(mhp_train, aes(x = rooms, y = price, colour = garage)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-habitaciones, garaje'),
  
  nrow =1
)

Representando boxplot por distrito se ve claramente que los barrios ricos como Salamanca o Chamberí, a parte de tener precios medios más altos, tienen mayor número de viviendas con precios atípicos (por encima), mientras que en el número de habitaciones no se aprecia una diferencia tan notoria.

ggplot(mhp_train, aes(district, price)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Precio", title = "Boxplot de precios por distrito") +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(district, rooms)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Habitaciones", title = "Boxplot de habitaciones por distrito") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(district, m2)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Metros cuadrados", title = "Boxplot de superficie por distrito") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

En las siguientes figuras se representa el precio medio por distrito y como varía con respecto a otras variables de interés.

mhp_train %>%
  group_by(district, m2) %>%
  summarize(avg_price = mean(price)) %>%
  ggplot(aes(x = m2, y = avg_price)) + 
  geom_point(size = 0.5) +
  facet_wrap(~ district) + 
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))

mhp_train %>%
  group_by(district, exterior) %>% 
  summarise(avg_price = mean(price)) %>%
  ggplot(aes(x=district, y=avg_price, fill=exterior)) +
  geom_bar(stat = "identity", position = "dodge") + 
  scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  ggtitle("Precio medio por distrito y exterior") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) + 
  labs(x = "Distrito", y = "Precio medio") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

mhp_train %>%
  group_by(district, elevator) %>% 
  summarise(avg_price = mean(price)) %>%
  ggplot(aes(x=district, y=avg_price, fill=elevator)) +
  geom_bar(stat = "identity", position = "dodge") + 
  scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  ggtitle("Precio medio por distrito y ascensor") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) + 
  labs(x = "Distrito", y = "Precio medio") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

mhp_train %>%
  group_by(district, garage) %>% 
  summarise(avg_price = mean(price)) %>%
  ggplot(aes(x=district, y=avg_price, fill=garage)) +
  geom_bar(stat = "identity", position = "dodge") + 
  scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  ggtitle("Precio medio por distrito y garaje") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) + 
  labs(x = "Distrito", y = "Precio medio") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

5 Procesado de variables cualitativas

mhp_train %>%
  group_by(house_type) %>%
  summarise(precio_medio = mean(price)) %>%
  arrange(precio_medio)
## # A tibble: 26 × 2
##    house_type  precio_medio
##    <fct>              <dbl>
##  1 sotano           279031.
##  2 planta -1        298188.
##  3 semi-sotano      302916.
##  4 bajo             316306.
##  5 planta 11        458019.
##  6 planta 1         540375.
##  7 planta 2         595355.
##  8 entreplanta      598834.
##  9 planta 10        605368.
## 10 planta 4         605711.
## # … with 16 more rows

Para el posterior análisis transformamos la columna house_type, agrupando en “piso” (independientemente de la planta en la que esté), “casa” (casa o chalet) y “sotano”.

mhp_train <- mhp_train %>%
  mutate(tipo_casa = case_when(
    grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
    grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
    TRUE ~ "piso"
  ))
mhp_train$tipo_casa <- as.factor(mhp_train$tipo_casa)

También se agruparán los distritos en función del precio medio del mismo: caro, medio o barato.

mhp_train <- mhp_train %>%
  group_by(district) %>%
  mutate(precio_distrito = mean(price)) %>%
  ungroup()

# Calcular los percentiles
percentiles <- quantile(mhp_train$precio_distrito, probs = c(0, 1/3, 2/3, 1))

# Asignar etiquetas
mhp_train$precio_distrito <- cut(mhp_train$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)

Además, de cara a desarrollar algunos modelos, se creará una variable target binaria que dividirá las casas entre “cara” y “barata” (1/0).

mhp_train$precio_bin <- ifelse(mhp_train$price > median(mhp_train$price), "cara", "barata")
mhp_train$precio_bin = ifelse(mhp_train$precio_bin == "cara", 1, 0)
mhp_train$precio_bin = factor(mhp_train$precio_bin, levels = c(1,0))
# Aplicamos todos los cambios también a test.

# Categorizamos house_type
mhp_test <- mhp_test %>%
  mutate(tipo_casa = case_when(
    grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
    grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
    TRUE ~ "piso"
  ))
mhp_test$tipo_casa <- as.factor(mhp_test$tipo_casa)

# Categorizamos los distritos
mhp_test <- mhp_test %>%
  group_by(district) %>%
  mutate(precio_distrito = mean(price)) %>%
  ungroup()

percentiles <- quantile(mhp_test$precio_distrito, probs = c(0, 1/3, 2/3, 1))

mhp_test$precio_distrito <- cut(mhp_test$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)

# Creamos variable target
mhp_test$precio_bin <- ifelse(mhp_test$price > median(mhp_test$price), "cara", "barata")
mhp_test$precio_bin = ifelse(mhp_test$precio_bin == "cara", 1, 0)
mhp_test$precio_bin = factor(mhp_test$precio_bin, levels = c(1,0))

6 Modelo de regresión lineal

Con el dataset ya preparado, se crea un modelo de regresión lineal.

En este caso, tras hacer diversas pruebas se ha optado por introducir sólo 3 variables (tipo_casa, m2 y precio_distrito), alcanzando un R2 de 0.7553.

Esta cifra podía incrementarse algo, sin embargo es a costa de introducir más variables y por tanto más complejidad al modelo, haciéndolo menos explicable.

lm_fit <- lm(price ~ m2+tipo_casa+precio_distrito, data=mhp_train)
summary(lm_fit)
## 
## Call:
## lm(formula = price ~ m2 + tipo_casa + precio_distrito, data = mhp_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3277290  -137875   -17495    95705  5483233 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -915949.86   30574.68 -29.958   <2e-16 ***
## m2                      6052.79      51.28 118.027   <2e-16 ***
## tipo_casapiso         673879.37   27865.36  24.183   <2e-16 ***
## tipo_casasotano       606283.79   28773.95  21.071   <2e-16 ***
## precio_distritomedio  101183.27   10337.15   9.788   <2e-16 ***
## precio_distritocaro   373747.30   11419.31  32.729   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 372600 on 7726 degrees of freedom
## Multiple R-squared:  0.7553, Adjusted R-squared:  0.7552 
## F-statistic:  4771 on 5 and 7726 DF,  p-value: < 2.2e-16
coef(lm_fit)
##          (Intercept)                   m2        tipo_casapiso 
##          -915949.864             6052.795           673879.371 
##      tipo_casasotano precio_distritomedio  precio_distritocaro 
##           606283.787           101183.273           373747.298
residuals=lm_fit$residuals
autoplot(lm_fit)

Observando los residuos se ve como hay cierta heterocedasticidad y en la gráfica QQ no se ajustan a la normal, por lo que el módelo no es bueno.

# Realizar predicciones en la partición test
predictions <- predict(lm_fit, newdata = mhp_test)

# Comparar las predicciones con los valores reales
rmse <- rmse(predictions, mhp_test$price)
mae <- mae(predictions, mhp_test$price)
r2 <- R2(predictions, mhp_test$price)

# Imprimir las métricas de evaluación
cat(paste0("RMSE: ", round(rmse, 2), "\n"))
## RMSE: 398263.79
cat(paste0("MAE: ", round(mae, 2), "\n"))
## MAE: 199633.64
cat(paste0("R-squared: ", round(r2, 2)))
## R-squared: 0.73

Un R2 de 0.73 (bastante similar al calculado en train) indica que se explica con esa probabilidad la variabilidad de la variable dependiente, por lo tanto podríamos considerarlo razonablemente útil, sin embargo sería preferible mejorar los valores del error de predicción promedio o el cuadrático.

7 Conclusiones preliminares

El análisis realizado hasta el momento ha servido para familiarizarse con el mercado inmobiliario de Madrid (o al menos parte de él) y conocer como se comportan algunas de las variables más relevantes en el valor de las viviendas.

Con estos conocimientos se abordarán nuevos modelos para intentar conseguir mejores predicciones y simultaneamente comprender y desarrollar estos mismos modelos.

str(mhp_train)
## tibble [7,732 × 12] (S3: tbl_df/tbl/data.frame)
##  $ price          : num [1:7732] 499000 385000 950000 460000 335000 580000 940000 85300 399000 125000 ...
##  $ house_type     : Factor w/ 26 levels "bajo","casa",..: 21 16 20 10 18 19 16 16 6 1 ...
##  $ exterior       : Factor w/ 2 levels "1","0": 1 1 1 1 1 1 1 1 1 1 ...
##  $ rooms          : num [1:7732] 4 1 3 3 4 3 2 3 3 2 ...
##  $ m2             : num [1:7732] 131 60 160 114 122 120 124 75 86 62 ...
##  $ elevator       : Factor w/ 2 levels "1","0": 1 2 1 1 1 1 1 1 2 2 ...
##  $ garage         : Factor w/ 2 levels "1","0": 2 2 2 1 2 2 2 2 2 2 ...
##  $ neighborhood   : Factor w/ 341 levels "12 de Octubre-Orcasur",..: 76 17 8 284 318 265 264 316 192 11 ...
##  $ district       : Factor w/ 21 levels "Arganzuela","barajas",..: 17 12 7 13 18 3 3 21 1 19 ...
##  $ tipo_casa      : Factor w/ 3 levels "casa","piso",..: 2 2 2 2 2 2 2 2 3 3 ...
##  $ precio_distrito: Factor w/ 3 levels "barato","medio",..: 2 3 3 1 1 3 3 1 1 1 ...
##  $ precio_bin     : Factor w/ 2 levels "1","0": 1 1 1 1 2 1 1 2 1 2 ...

8 Aprendizaje no supervisado

8.1 K-Means

El algoritmo K-Means es un método de aprendizaje automático no supervisado utilizado para el análisis de agrupamiento o “clustering”. El objetivo principal de este algoritmo es dividir los datos en K grupos (clusters), de modo que los puntos dentro de un mismo grupo sean lo más similares posible.

Por como funciona este método, es habitual utilizarlo en la etapa de EDA para agrupar y etiquetar observaciones con un criterior más científico que el basado en el conocimiento personal del dominio, más que para hacer predicciones (aunque también puede hacerse).

En nuestro lo utilizaremos para comprobar si hay alguna forma mejor de agrupar las casas de como lo hemos hecho anteriormente.

mhp_train_kmeans <- subset(mhp_train, select = c("m2", "rooms", "elevator", "garage", "exterior"))

mhp_train_kmeans$exterior <- as.numeric(mhp_train_kmeans$exterior)
mhp_train_kmeans$elevator <- as.numeric(mhp_train_kmeans$elevator)
mhp_train_kmeans$garage <- as.numeric(mhp_train_kmeans$garage)

# Normalizamos los datos
mhp_train_kmeans <- scale(mhp_train_kmeans)

Tras normalizar los datos, aplicamos k-means sólo para las variables numéricas y mediante los siguientes test decidiremos cual es el número óptimo de clusters.

fviz_nbclust(mhp_train_kmeans, kmeans, method = "wss")

fviz_nbclust(mhp_train_kmeans, kmeans, method = "silhouette")

Utilizando la regla del codo en el gráfico wss podríamos decir que está en el 6, algo que confirma el segundo gráfico. Por lo tanto, aplicamos el modelo con K = 6

# Definir el número de grupos
k <- 6

# Ejecutar k-means
grupos <- kmeans(mhp_train_kmeans, k)

mhp_train$cluster <- as.factor(grupos$cluster)

# Graficar la relación entre las características y los grupos
ggplot(mhp_train, aes(x=m2, y=price, color=cluster)) + 
  geom_point() +
  ggtitle("Grupos de casas en venta") +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Tamaño en m2") +
  ylab("Precio")

fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "euclid",repel = TRUE,star.plot = TRUE)

fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "norm")

Graficamente es difícil diferenciar las agrupaciones, por lo que en la tabla inferior obtenemos los precios medios de cada cluster, donde se puede apreciar que la gran mayoría de agrupaciones son coherentes (salvo tal vez el 2 y el 3 que son muy similares).

También se muestra el tipo de casa y el distrito predominante en cada cluster.

# Precio medio por cluster
mhp_train %>%
  group_by(cluster) %>%
  summarise(precio_medio = mean(price))
## # A tibble: 6 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            184203.
## 2 2            640443.
## 3 3            617779.
## 4 4            223476.
## 5 5           2228444.
## 6 6            355840.
# Tipo de casa predominante por cluster
mhp_train %>%
  group_by(cluster, tipo_casa) %>%
  summarise(frecuencia = n()) %>%
  arrange(cluster, desc(frecuencia)) %>%
  slice(1) %>%
  select(cluster, tipo_casa_predominante = tipo_casa)
## # A tibble: 6 × 2
## # Groups:   cluster [6]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       sotano                
## 2 2       piso                  
## 3 3       piso                  
## 4 4       piso                  
## 5 5       piso                  
## 6 6       piso
# Distrito predominante por cluster
mhp_train %>%
  group_by(cluster, district) %>%
  summarise(frecuencia = n()) %>%
  arrange(cluster, desc(frecuencia)) %>%
  slice(1) %>%
  select(cluster, distrito_predominante = district)
## # A tibble: 6 × 2
## # Groups:   cluster [6]
##   cluster distrito_predominante
##   <fct>   <fct>                
## 1 1       tetuan               
## 2 2       hortaleza            
## 3 3       barrio de salamanca  
## 4 4       carabanchel          
## 5 5       barrio de salamanca  
## 6 6       barrio de salamanca

Sabemos que tanto la regla del codo como cualquier otro método para definir el número de clusters no siempre es adecuada, por lo que a continuación se prueba con un bucle que prueba con número de clusters menores que 6 y muestra las agrupaciones por precio y tipo de piso anteriores.

# Ejecuta k-means para diferentes valores de K
for (k in 2:6) {
  # Ejecutar k-means
  grupos <- kmeans(mhp_train_kmeans, centers = k)
  mhp_train$cluster <- as.factor(grupos$cluster)
  
  print(k)
  
  print(mhp_train %>%
  group_by(cluster) %>%
  summarise(precio_medio = mean(price)))
  
  print(mhp_train %>%
  group_by(cluster, tipo_casa) %>%
  summarise(frecuencia = n()) %>%
  arrange(cluster, desc(frecuencia)) %>%
  slice(1) %>%
  select(cluster, tipo_casa_predominante = tipo_casa))
}
## [1] 2
## # A tibble: 2 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            394442.
## 2 2           1094473.
## # A tibble: 2 × 2
## # Groups:   cluster [2]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## [1] 3
## # A tibble: 3 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            402986.
## 2 2            605000.
## 3 3           1941976.
## # A tibble: 3 × 2
## # Groups:   cluster [3]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       piso                  
## [1] 4
## # A tibble: 4 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            579636.
## 2 2            352563.
## 3 3            665255.
## 4 4            189057.
## # A tibble: 4 × 2
## # Groups:   cluster [4]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       piso                  
## 4 4       sotano                
## [1] 5
## # A tibble: 5 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            435798.
## 2 2            241929.
## 3 3            189057.
## 4 4            665853.
## 5 5           1005091.
## # A tibble: 5 × 2
## # Groups:   cluster [5]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       sotano                
## 4 4       piso                  
## 5 5       piso                  
## [1] 6
## # A tibble: 6 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1           1179261.
## 2 2            652498.
## 3 3            355840.
## 4 4            217255.
## 5 5           2356279.
## 6 6            451222.
## # A tibble: 6 × 2
## # Groups:   cluster [6]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       piso                  
## 4 4       piso                  
## 5 5       piso                  
## 6 6       piso

Observando las diferencias entre los diferentes modelos, el número de clusters que mejor diferencia las casas por su precio es 3, algo que encaja con la presunción que inicial que hicimos para dividir los tipos de casa en 3 tipos.

Si observamos además el tipo de casas más común por cluster, casi siempre es piso (salvo en los casos de clusters con muy poco precio medio, que es sótano), lo cual aunque tiene sentido porque este tipo es claramente mayoritario en todo el dataset, no nos permite apreciar a simple vista donde se acumulan el resto de viviendas.

Por lo tanto, entre los datos extraídos del análisis de clusters y el conocimiento previo que teníamos del dominio, tiene sentido mantener los 3 tipos de viviendas principales como variable para los siguientes modelos (sótano, piso y casa).

9 Técnicas de reducción de la dimensionalidad

9.1 PCA

El Análisis de Componentes Principales es una técnica de reducción de la dimensionalidad que se utiliza para transformar un conjunto de datos de múltiples variables en un nuevo conjunto de datos con un número menor de variables cuyo principal objetivo es simplificar la representación de datos mientras se conserva la mayor cantidad de información y variabilidad posible identificando patrones y estructuras subyacentes, reduciendo la multicolinealidad y el sobreajuste.

PAra nuestro dataset, seleccionaremos las columnas numéricas, descartando la columna precio puesto que será de alguna manera nuestra variable dependiente. Además creamos 3 variables dummy usando variables categóricas de tipo 0|1.

train_pca <- mhp_train[,c(3:7,10,11,12)]
dummy_garaje <- model.matrix(~0 + garage,data = train_pca)
dummy_exterior <- model.matrix(~0 + exterior, data = train_pca)
dummy_ascensor <- model.matrix(~0 + elevator, data = train_pca)
dummy_tipo_casa <- model.matrix(~0 + tipo_casa, data = train_pca)
dummy_distrito <- model.matrix(~0 + precio_distrito, data = train_pca)
dummy_data <- cbind(train_pca[, c("rooms", "m2")],dummy_exterior,dummy_garaje,dummy_ascensor,dummy_distrito,dummy_tipo_casa)
prcomp(dummy_data)
## Standard deviations (1, .., p=14):
##  [1] 1.025015e+02 9.457105e-01 7.267966e-01 6.257613e-01 6.127629e-01
##  [6] 5.782676e-01 4.571490e-01 4.037562e-01 1.716715e-01 2.268386e-15
## [11] 1.853843e-15 1.003401e-15 8.037840e-16 4.438592e-16
## 
## Rotation (n x k) = (14 x 14):
##                                 PC1           PC2          PC3          PC4
## rooms                 -8.797027e-03 -0.9916184565 -0.014616928  0.094846395
## m2                    -9.999557e-01  0.0085416356  0.001971791 -0.001450128
## exterior1             -5.771447e-04 -0.0453751425 -0.009449951 -0.063478191
## exterior0              5.771447e-04  0.0453751425  0.009449951  0.063478191
## garage1               -1.538147e-03  0.0421269404 -0.223561300  0.037355645
## garage0                1.538147e-03 -0.0421269404  0.223561300 -0.037355645
## elevator1             -4.274374e-04  0.0241411854 -0.476833313  0.342257253
## elevator0              4.274374e-04 -0.0241411854  0.476833313 -0.342257253
## precio_distritobarato  1.400700e-03 -0.0385563529  0.200000119 -0.359200946
## precio_distritomedio   6.834450e-05 -0.0051663464 -0.051807433  0.201186914
## precio_distritocaro   -1.469044e-03  0.0437226993 -0.148192687  0.158014031
## tipo_casacasa         -8.085078e-04  0.0005892036  0.067961955 -0.037840260
## tipo_casapiso          2.818685e-05 -0.0469938502 -0.466029098 -0.505128722
## tipo_casasotano        7.803209e-04  0.0464046466  0.398067143  0.542968982
##                                 PC5           PC6          PC7          PC8
## rooms                 -0.0132115607 -0.0291253564  0.074126853 -0.029792435
## m2                     0.0009596296  0.0005665667 -0.002313108 -0.001386419
## exterior1             -0.1552630611 -0.0333936354 -0.327138796  0.601052285
## exterior0              0.1552630611  0.0333936354  0.327138796 -0.601052285
## garage1               -0.5627134868 -0.1309300044  0.333444368  0.033280690
## garage0                0.5627134868  0.1309300044 -0.333444368 -0.033280690
## elevator1              0.0482119831 -0.1533456743 -0.318624110 -0.141577832
## elevator0             -0.0482119831  0.1533456743  0.318624110  0.141577832
## precio_distritobarato -0.2366848900 -0.4502445095 -0.365221920 -0.322493702
## precio_distritomedio  -0.1969164022  0.7585991702 -0.088292084 -0.036423761
## precio_distritocaro    0.4336012922 -0.3083546607  0.453514004  0.358917463
## tipo_casacasa         -0.0381883034  0.0229829599  0.061657352 -0.001868362
## tipo_casapiso          0.1388358808  0.1326235059  0.004437444 -0.006913037
## tipo_casasotano       -0.1006475774 -0.1556064657 -0.066094796  0.008781399
##                                PC9          PC10          PC11          PC12
## rooms                  0.000862928  5.877890e-16  1.174221e-15  7.345674e-17
## m2                    -0.001036478 -4.349751e-20  3.574359e-18 -4.762673e-18
## exterior1              0.017792902  8.077407e-03  2.416838e-02 -9.348594e-03
## exterior0             -0.017792902  8.077407e-03  2.416838e-02 -9.348594e-03
## garage1               -0.027676182  6.464949e-03  5.102210e-03 -8.893342e-02
## garage0                0.027676182  6.464949e-03  5.102210e-03 -8.893342e-02
## elevator1              0.086614598  6.353245e-01 -3.028358e-01 -6.708710e-02
## elevator0             -0.086614598  6.353245e-01 -3.028358e-01 -6.708710e-02
## precio_distritobarato -0.004863213 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritomedio  -0.010432322 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritocaro    0.015295535 -9.943406e-02 -8.213108e-02 -5.584679e-01
## tipo_casacasa          0.809211511 -2.329879e-01 -5.148217e-01  1.145226e-01
## tipo_casapiso         -0.393939361 -2.329879e-01 -5.148217e-01  1.145226e-01
## tipo_casasotano       -0.415272150 -2.329879e-01 -5.148217e-01  1.145226e-01
##                                PC13          PC14
## rooms                 -9.651041e-17 -5.224615e-18
## m2                     4.139246e-18 -4.009145e-19
## exterior1              5.278108e-02 -7.046115e-01
## exterior0              5.278108e-02 -7.046115e-01
## garage1               -6.995888e-01 -5.097578e-02
## garage0               -6.995888e-01 -5.097578e-02
## elevator1              1.228500e-02 -1.293890e-03
## elevator0              1.228500e-02 -1.293890e-03
## precio_distritobarato  6.884856e-02  8.609917e-03
## precio_distritomedio   6.884856e-02  8.609917e-03
## precio_distritocaro    6.884856e-02  8.609917e-03
## tipo_casacasa         -1.877163e-02 -2.325502e-02
## tipo_casapiso         -1.877163e-02 -2.325502e-02
## tipo_casasotano       -1.877163e-02 -2.325502e-02

Podemos observar que en la PC1 lleva un peso casi de 100 la variable m2, mientras que en PC2 pasa lo mismo con la variable rooms. En el resto de PCs el análisis de componentes principales está algo más repartido.

summary(prcomp(dummy_data))
## Importance of components:
##                             PC1     PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     102.5015 0.94571 0.72680 0.62576 0.61276 0.57827 0.45715
## Proportion of Variance   0.9997 0.00009 0.00005 0.00004 0.00004 0.00003 0.00002
## Cumulative Proportion    0.9997 0.99981 0.99986 0.99989 0.99993 0.99996 0.99998
##                            PC8    PC9      PC10      PC11      PC12      PC13
## Standard deviation     0.40376 0.1717 2.268e-15 1.854e-15 1.003e-15 8.038e-16
## Proportion of Variance 0.00002 0.0000 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  1.00000 1.0000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
##                             PC14
## Standard deviation     4.439e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

Al hacer un summary de las variables vemos que la proporción de varianza explicada por la primera componente es de un 99%, es decir, que es la única relevante. Si nos fijamos vemos que la PC1 está explicada casi al 100% por los m2. En nuestra tabla vemos que los m2 tienen valores mucho más altos que el resto de variables. Por lo tanto, lo que hacemos a continuación es un análisis pero escalando las variables.

prcomp(dummy_data, scale = T)
## Standard deviations (1, .., p=14):
##  [1] 1.788828e+00 1.554767e+00 1.391248e+00 1.320487e+00 1.255039e+00
##  [6] 1.218340e+00 9.006738e-01 7.679348e-01 4.930742e-01 2.297198e-14
## [11] 4.169026e-15 2.788355e-15 1.325148e-15 8.186691e-16
## 
## Rotation (n x k) = (14 x 14):
##                               PC1          PC2         PC3         PC4
## rooms                 -0.32118124  0.185174155  0.19374853 -0.18231760
## m2                    -0.39188093  0.178866531  0.28992752 -0.15926158
## exterior1             -0.22525556  0.323113180 -0.41045579  0.16000993
## exterior0              0.22525556 -0.323113180  0.41045579 -0.16000993
## garage1               -0.39123327  0.104563202  0.02079773  0.20853492
## garage0                0.39123327 -0.104563202 -0.02079773 -0.20853492
## elevator1             -0.28418360 -0.450986811 -0.07378472  0.27110372
## elevator0              0.28418360  0.450986811  0.07378472 -0.27110372
## precio_distritobarato  0.20341246  0.208545636 -0.33198104  0.01261772
## precio_distritomedio  -0.02774634 -0.002015776  0.03875669  0.15703618
## precio_distritocaro   -0.18188148 -0.213762808  0.30358073 -0.17520956
## tipo_casacasa         -0.14205765  0.358486679  0.29646454 -0.20249803
## tipo_casapiso         -0.16858177 -0.251362769 -0.39199358 -0.48519863
## tipo_casasotano        0.22431380  0.124024528  0.29030539  0.56945631
##                               PC5          PC6          PC7          PC8
## rooms                  0.17534372 -0.091653082  0.472849806  0.496683083
## m2                     0.12920863 -0.038990789  0.209138336  0.106709918
## exterior1              0.31433452 -0.162497819 -0.145528928 -0.064746601
## exterior0             -0.31433452  0.162497819  0.145528928  0.064746601
## garage1               -0.40131726  0.276379019 -0.203710793  0.095975695
## garage0                0.40131726 -0.276379019  0.203710793 -0.095975695
## elevator1              0.11627307  0.012675714  0.278402861 -0.212854201
## elevator0             -0.11627307 -0.012675714 -0.278402861  0.212854201
## precio_distritobarato -0.04064647  0.512859981  0.434958282 -0.014241486
## precio_distritomedio  -0.38476634 -0.695032433  0.035559829  0.020880654
## precio_distritocaro    0.43936573  0.186858173 -0.486898651 -0.006820754
## tipo_casacasa         -0.08089606 -0.001818982  0.145439408 -0.763400011
## tipo_casapiso         -0.15088701 -0.039033227 -0.059791130  0.108495318
## tipo_casasotano        0.18370301  0.040471673  0.007300261  0.171061449
##                               PC9          PC10          PC11          PC12
## rooms                 -0.53136448  1.871326e-14  6.573067e-16 -8.713797e-16
## m2                     0.79477851  5.562477e-15  5.135097e-16  5.962879e-16
## exterior1             -0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## exterior0              0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## garage1               -0.06271669 -7.002484e-04 -1.864879e-03  2.225429e-02
## garage0                0.06271669 -7.002484e-04 -1.864879e-03  2.225429e-02
## elevator1             -0.02034701 -1.325824e-03 -1.363473e-02  3.145845e-02
## elevator0              0.02034701 -1.325824e-03 -1.363473e-02  3.145845e-02
## precio_distritobarato  0.09621994 -7.055211e-02  4.768696e-02 -5.771886e-01
## precio_distritomedio   0.02039646 -7.038638e-02  4.757494e-02 -5.758327e-01
## precio_distritocaro   -0.12064810 -6.816645e-02  4.607446e-02 -5.576714e-01
## tipo_casacasa         -0.22215925 -2.483203e-01 -1.912383e-03  3.007261e-02
## tipo_casapiso          0.04299775 -6.861954e-01 -5.284580e-03  8.310109e-02
## tipo_casasotano        0.03813232 -6.729672e-01 -5.182706e-03  8.149911e-02
##                                PC13          PC14
## rooms                 -4.757936e-17  3.988302e-17
## m2                    -1.440589e-18 -3.087024e-16
## exterior1              1.099278e-02  3.536496e-04
## exterior0              1.099278e-02  3.536496e-04
## garage1                2.157713e-02 -7.064242e-01
## garage0                2.157713e-02 -7.064242e-01
## elevator1             -7.059752e-01 -2.053508e-02
## elevator0             -7.059752e-01 -2.053508e-02
## precio_distritobarato -2.595459e-02 -1.903173e-02
## precio_distritomedio  -2.589362e-02 -1.898702e-02
## precio_distritocaro   -2.507696e-02 -1.838819e-02
## tipo_casacasa          1.806856e-03  1.253756e-03
## tipo_casapiso          4.992972e-03  3.464565e-03
## tipo_casasotano        4.896719e-03  3.397776e-03

Tras escalar las variables podemos observar que están mucho más repartidos los pesos en las componentes principales debido a que estás asegurando que todas las variables tengan la misma importancia en el análisis, independientemente de sus escalas originales y unidades de medida. Con esto evitas el dominio de variables con mayor magnitud y mejoras la interpretación de resultados.

pca_result <- prcomp(dummy_data, center = TRUE, scale = TRUE)
summary(pca_result)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5    PC6     PC7
## Standard deviation     1.7888 1.5548 1.3912 1.3205 1.2550 1.2183 0.90067
## Proportion of Variance 0.2286 0.1727 0.1383 0.1245 0.1125 0.1060 0.05794
## Cumulative Proportion  0.2286 0.4012 0.5395 0.6640 0.7765 0.8826 0.94051
##                            PC8     PC9      PC10      PC11      PC12      PC13
## Standard deviation     0.76793 0.49307 2.297e-14 4.169e-15 2.788e-15 1.325e-15
## Proportion of Variance 0.04212 0.01737 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  0.98263 1.00000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
##                             PC14
## Standard deviation     8.187e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

Observamos que con PC6 ya llegamos al 88% de la proporción de varianza explicada pero vemos algunos saltos interesantes. Para ver cuantas componentes explican el modelo hacemos la regla de codo, método heurístico para determinar el número óptimo de clusters en un conjunto de datos.

Al observar el gráfico vemos que la forma del codo sería en la 6|7 aunque el gráfico no es muy claro.

summary_pca <- summary(pca_result)
var_exp <- summary_pca$importance[2,]

barplot(var_exp, col = "#0BB363",
        main = "Gráfico PCA - Variance explained",
        xlab = "Principal Components",
        ylab = "Proportion of Variance Explained")

lines(var_exp, col = "#DC143C", type = "l", lwd = 2)

points(var_exp, col = "#DC143C", pch = 20, cex = 1.2)

A modo de ejercicio para comprender el funcionamiento de PCA puede resultar útil hacer estos cálculos, sin embargo para el dataset en cuestión ya se pueden formar modelos robustos con muy pocas variables (los m2 de superficie, por ejemplo, ya son suficientes para predecir un gran porcentaje de precios).

10 Aprendizaje supervisado

10.1 GLM

En aprendizaje automático, GLM (Generalized Linear Models) se refiere a una clase de modelos lineales que se utilizan para predecir una variable respuesta a partir de una o más variables explicativas. Tiene tres componentes: Variable respuesta, función de enlace y función lineal sistemática, y son muy útiles para utilizarse en clasificación, regresión y problemas de predicción.

Aplicado a nuestro dataset, lo utilizaremos para predecir el precio binario de las casas (cara/barata).

train_glm <- mhp_train[,c(3:7,10,11,12)]
test_glm <- mhp_test[,c(3:7,10,11,12)]
train_glm[,2:3]<- scale(train_glm[,2:3])
test_glm[,2:3]<- scale(test_glm[,2:3])

Tras seleccionar las variables y normalizarlas para que todas tengan el mismo peso, desarrollamos el modelo con la variable de respuesta “precio_bin” y las variables independientes “exterior”, “garage”, “elevator”, “rooms”, “precio_distrito” y “tipo_casa”, usando una distribución binomial.

glm_mhp = glm(formula = precio_bin ~ exterior + garage + elevator + rooms + m2 + precio_distrito + tipo_casa,
                 data = train_glm, 
                 family = binomial)
summary(glm_mhp)
## 
## Call:
## glm(formula = precio_bin ~ exterior + garage + elevator + rooms + 
##     m2 + precio_distrito + tipo_casa, family = binomial, data = train_glm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5396  -0.1705   0.0219   0.2836   5.0687  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.99261    0.62333  -1.592  0.11129    
## exterior0            -0.35269    0.13261  -2.660  0.00782 ** 
## garage0               0.79113    0.10737   7.368 1.73e-13 ***
## elevator0             1.67266    0.12813  13.054  < 2e-16 ***
## rooms                -0.01548    0.08151  -0.190  0.84940    
## m2                   -6.20689    0.23649 -26.246  < 2e-16 ***
## precio_distritomedio -2.52394    0.11266 -22.403  < 2e-16 ***
## precio_distritocaro  -4.33905    0.14476 -29.974  < 2e-16 ***
## tipo_casapiso         0.34751    0.62328   0.558  0.57716    
## tipo_casasotano       1.01944    0.62491   1.631  0.10282    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10718.8  on 7731  degrees of freedom
## Residual deviance:  3535.9  on 7722  degrees of freedom
## AIC: 3555.9
## 
## Number of Fisher Scoring iterations: 8

En este caso, todos los coeficientes tienen valores “p” muy pequeños, lo que indica que todas las variables independientes tienen una relación significativa con la variable dependiente “precio_bin” en el modelo de regresión logística ajustado.

Siendo el “Intercept” negativo, se puede apreciar como el precio depende de forma proporcional con el número de habitaciones, el tipo de piso y los metros cuadrados (sobre todo). Todo esperable, así como la dependencia inversa con aquellas viviendas que no tienen ascensor o graje.

El modelo explica una cantidad significativa de la varianza en la variable respuesta, con un AIC de 4818.

head(predict(glm_mhp))
##          1          2          3          4          5          6 
## -2.6264925  1.5662145 -6.1858748  0.1476553  0.4424557 -3.7636044

Esto sirve para predecir la probabilidad logit ajustado previamente y para mostrar las primeras seis predicciones. Cada número representa la probabilidad logit de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.

ajustados <- fitted(glm_mhp)
head(fitted(glm_mhp))
##           1           2           3           4           5           6 
## 0.067452746 0.827243283 0.002054072 0.536846908 0.608844012 0.022673932

Aquí obtenemos las probabilidades ajustadas que representan las probabilidades estimadas de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.

Gráficos de residuales: Estos gráficos permiten evaluar la calidad del ajuste del modelo. Puedes trazar residuales de Pearson, deviance, o residuales estandarizados frente a las predicciones ajustadas o valores ajustados para examinar si hay patrones en los residuales. Un buen ajuste del modelo no mostrará patrones claros en los residuales.

residuales <- residuals(glm_mhp, type = "pearson")
head(residuales)
##           1           2           3           4           5           6 
## -0.26894557 -2.18826117 -0.04536849 -1.07662124  0.80153404 -0.15231535
data_residuales <- data.frame(Residuales = residuales, Ajustados = ajustados)
ggplot(data_residuales, aes(x = Ajustados, y = Residuales)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Gráfico de residuales", x = "Valores ajustados", y = "Residuales") + ylim(-10,10)

predicciones_glm <- ifelse(test = glm_mhp$fitted.values > 0.5, no = 0, yes = 1)
matriz_confusion_glm <- table(glm_mhp$model$precio_bin, predicciones_glm,
                          dnn = c("observaciones", "predicciones"))
matriz_confusion_glm
##              predicciones
## observaciones    0    1
##             1 3449  414
##             0  329 3540

Finalmente generamos predicciones binarias basadas en las probabilidades ajustadas del modelo lineal generalizado y creamos una matriz de confusión para evaluar el rendimiento del modelo donde se muestra la distribución de las predicciones correctas e incorrectas en función de las observaciones reales y las predichas.

La matriz se lee de la siguiente manera:

  • Verdaderos negativos (VN): La esquina superior izquierda (3265) representa el número de observaciones clasificadas correctamente como “0” (clase negativa). Esto significa que el modelo predijo correctamente 3265 casos como “0”.

  • Falsos positivos (FP): La esquina superior derecha (598) representa el número de observaciones que en realidad eran “0” pero fueron clasificadas incorrectamente como “1” (clase positiva) por el modelo. El modelo cometió un error en 598 casos al clasificarlos como “1”.

  • Falsos negativos (FN): La esquina inferior izquierda (501) representa el número de observaciones que en realidad eran “1” pero fueron clasificadas incorrectamente como “0” por el modelo. El modelo cometió un error en 501 casos al clasificarlos como “0”.

  • Verdaderos positivos (VP): La esquina inferior derecha (3368) representa el número de observaciones clasificadas correctamente como “1”. Esto significa que el modelo predijo correctamente 3368 casos como “1”.

10.2 KNN

El modelo KNN (K-Nearest Neighbors) es un algoritmo de aprendizaje supervisado que se utiliza tanto para clasificación como para regresión. Es un método basado en instancias que no hace supuestos sobre la forma funcional de la relación entre las variables independientes y la variable dependiente. KNN utiliza la información de los vecinos más cercanos para realizar predicciones. Para ello calcula distancias como la Euclidiana, Manhattan, Minkowski, etc…

La idea detrás de KNN es simple, una observación se clasifica o se le asigna un valor basado en las características de sus K vecinos más cercanos.

Algunas ventajas del algoritmo KNN son su simplicidad o su capacidad para manejar relaciones no lineales entre variables. Sin embargo también tiene algunas desventajas, como su sensibilidad a la maldición de la dimensionalidad, el efecto de las variables irrelevantes y el costo computacional asociado con el cálculo de distancias en conjuntos de datos grandes.

Procedemos al cálculo de este modelo e intentamos ver si hay un número correcto de k-vecinos que podemos usar.

long = 15
accuracy = rep(0,long)
f1score = rep(0,long)
recall = rep(0,long)
precision = rep(0,long)
for (i in 1:long)
{
  prediccion_knn_cv =knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], 
                            k=i, cl=mhp_train$precio_bin)
  accuracy[i] = sum(prediccion_knn_cv == mhp_train$precio_bin) /nrow(mhp_train)
  recall[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & mhp_train$precio_bin == TRUE) / sum(mhp_train$precio_bin == TRUE)
  precision[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & prediccion_knn_cv == TRUE) / sum(prediccion_knn_cv == TRUE)
  f1score[i] = 2*precision[i]*recall[i]/(precision[i]+recall[i])
}
resultados_knn = as.data.frame(cbind(accuracy,f1score,precision,recall))
resultados_knn = resultados_knn %>% mutate(index=as.factor(seq(1:long)))

max(resultados_knn$f1score)
## [1] NaN
which.max(resultados_knn$f1score)
## integer(0)
ggplot(data=resultados_knn,aes(x=index,y=accuracy)) + 
  geom_col(colour="cyan4",fill="cyan3")+
  ggtitle("Accuracy")

ggplot(data=resultados_knn,aes(x=index,y=f1score)) + 
  geom_col(colour="orange4",fill="orange3") +
  ggtitle("F1_score values")

Tratamos de encontrar el valor óptimo de k utilizando validación cruzada. Se calculan varias métricas de evaluación, como precisión, recall y F1-score, para cada valor de k y finalmente tratamos de graficar la precisión en función del valor de k, pero analizando el gráfico no observamos ninguna diferencia, así que cogemos 5, que es el tamaño por convención que se suele coger.

Pasamos a realizar predicciones en el conjunto de datos de entrenamiento y en el de test además de realizar la matriz de confusión.

# En train
prediccion_knn5_train =knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], 
                              k=5, cl=mhp_train$precio_bin)
confusionMatrix(table(prediccion_knn5_train,mhp_train$precio_bin), positive= "1")
## Confusion Matrix and Statistics
## 
##                      
## prediccion_knn5_train    1    0
##                     1 3233  531
##                     0  630 3338
##                                           
##                Accuracy : 0.8498          
##                  95% CI : (0.8417, 0.8577)
##     No Information Rate : 0.5004          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6997          
##                                           
##  Mcnemar's Test P-Value : 0.004026        
##                                           
##             Sensitivity : 0.8369          
##             Specificity : 0.8628          
##          Pos Pred Value : 0.8589          
##          Neg Pred Value : 0.8412          
##              Prevalence : 0.4996          
##          Detection Rate : 0.4181          
##    Detection Prevalence : 0.4868          
##       Balanced Accuracy : 0.8498          
##                                           
##        'Positive' Class : 1               
## 
#En test
prediccion_knn5_test=knn(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], mhp_test[,c("exterior","rooms","m2","elevator", "garage")],
                         k=5, cl=mhp_train$precio_bin)
confusionMatrix(table(prediccion_knn5_test,mhp_test$precio_bin), positive= "1")
## Confusion Matrix and Statistics
## 
##                     
## prediccion_knn5_test    1    0
##                    1 1577  250
##                    0  356 1683
##                                           
##                Accuracy : 0.8432          
##                  95% CI : (0.8314, 0.8546)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6865          
##                                           
##  Mcnemar's Test P-Value : 1.996e-05       
##                                           
##             Sensitivity : 0.8158          
##             Specificity : 0.8707          
##          Pos Pred Value : 0.8632          
##          Neg Pred Value : 0.8254          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4079          
##    Detection Prevalence : 0.4726          
##       Balanced Accuracy : 0.8432          
##                                           
##        'Positive' Class : 1               
## 

La matriz de confusión en train nos muestra la cantidad de predicciones correctas e incorrectas de cada clase (1 y 0) en el conjunto de entrenamiento.

  • Accuracy (exactitud): 0.8482. Es la proporción de predicciones correctas en el conjunto de entrenamiento. En otras palabras, el modelo KNN clasifica correctamente el 84.82% de las observaciones del conjunto de entrenamiento.

  • Kappa: 0.6963. El índice kappa mide la concordancia entre las predicciones y los valores reales. Un kappa de 1 indica una concordancia perfecta, mientras que un kappa de 0 indica una concordancia totalmente aleatoria. Nuestro kappa de 0.6963 muestra una buena concordancia entre las predicciones y los valores reales en el conjunto de entrenamiento.

  • Mcnemar’s Test P-Value: 0.003201. Esta prueba evalúa si hay una diferencia significativa entre las predicciones falsas positivas y las predicciones falsas negativas.

La matriz de confusión en test nos muestra valores muy parecidos en los indicadores anteriores lo que nos lleva a pensar que el modelo tiene un rendimiento consistente y generaliza bien lo datos ya que el modelo clasifica correctamente las observaciones y no parece sobreajustado.

Hemos realizado diferentes pruebas quitando variables y se observa que si quitamos la variable m2 nos dice que hay muchos empates. Sin embargo, quitando las otras variables el acierto en train y test varía.

Obteniendo más o menos un 84% de acierto creemos que nos puede servir para clasificar las nuevas casas que entrasen en el dataset, puesto que no hemos sido capaces de mejorar ese %.

10.3 Decision Trees

El modelo de Árbol de Decisión es un algoritmo de aprendizaje supervisado que funciona dividiendo iterativamente el conjunto de datos en subconjuntos basados en las características del mismo. Al final de este proceso, se genera un árbol con nodos y hojas que representan las divisiones y las decisiones.

Algunas características importantes de los árboles de decisión son la interpretabilidad, capacidad para manejar datos no lineales, el poco preprocesamiento de los datos, la sensibilidad a datos desequilibrados y el posible sobreajuste en el caso de hacer un arbol con demasiadas ramificaciones.

A continuación hacemos el modelo de arbol de decisión.

arbol <- rpart(precio_bin ~ garage + elevator + rooms + tipo_casa + precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")

Para evaluar el modelo utilizamos la función predict sobre el conjunto test usando el argumento type = “class” que indica que se espera una salida categórica creandose una tabla de contingencias “tab1” con las predicciones que coinciden o no respecto a las observaciones reales.

Calculamos y obtenemos el número 0.8396275 que representa la tasa de acierto (accuracy) del modelo. El Árbol de Decisión hizo predicciones correctas en casi el 84% de los casos en el conjunto de datos de prueba. Nos parece un porcentaje bastante aceptable.

tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1
## [1] 0.8396275

Siguiendo con la validación del modelo podemos usar las funciones printcp y plotcp que se utilizan para analizar la complejidad y el rendimiento usando el “pruning” o poda de ramas del arbol que a menudo reduce su complejidad y permite una mejor generalización del mismo disminuyendo el subreajuste.

La función printcp es una función que muestra una tabla que contiene información sobre la complejidad y la tasa de error en diferentes niveles de tamaño del árbol proporcionando una visión general de cómo cambia el rendimiento a medida que aumenta la complejidad. El valor de cp ha de ser tal que la tasa de error de validación cruzada sea mínima.

La función plotcp crea un gráfico de la tasa de error. Este gráfico es útil para visualizar cómo se comporta el rendimiento del árbol a medida que aumenta su tamaño y ayuda a decidir qué tamaño de árbol es el más apropiado para el problema en cuestión.

printcp(arbol)
## 
## Classification tree:
## rpart(formula = precio_bin ~ garage + elevator + rooms + tipo_casa + 
##     precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
## 
## Variables actually used in tree construction:
## [1] elevator        precio_distrito rooms          
## 
## Root node error: 3863/7732 = 0.49961
## 
## n= 7732 
## 
##         CP nsplit rel error  xerror      xstd
## 1 0.516438      0   1.00000 1.03055 0.0113762
## 2 0.044784      1   0.48356 0.48356 0.0097435
## 3 0.041677      3   0.39399 0.40098 0.0091108
## 4 0.016179      4   0.35232 0.35232 0.0086689
## 5 0.010000      6   0.31996 0.31996 0.0083418

Obtenidos los resultados, pasamos a analizarlos:

  • CP (Cost Complexity): Mide la complejidad asociado con cada nivel de poda en el árbol de decisión. Un CP más pequeño indica una complejidad menor.
  • nsplit: Número de divisiones (splits) realizadas hasta el momento en el árbol.
  • rel error: Error relativo en cada etapa de poda del árbol.
  • xerror: Error de validación cruzada (cross-validated) estimado para el árbol en cada etapa de poda.
  • xstd: Desviación estándar del error de validación cruzada (xerror).

En base a los resultados presentados, se puede observar que el árbol sin poda tiene un error relativo de 1.0 y un error de validación cruzada de 1.0189. A medida que se realizan más divisiones en el árbol el CP disminuye y también lo hace el error relativo y el error de validación cruzada.

Para seleccionar el árbol óptimo, generalmente se busca el CP más pequeño que tenga un error de validación cruzada dentro de una desviación estándar del mínimo error de validación cruzada.

arbol$cptable[which.min(arbol$cptable[, "xerror"]),
"CP"]
## [1] 0.01

El propósito es identificar el mejor árbol de decisión con respecto a su capacidad de generalización en datos no vistos. La función which.min() devuelve el índice del valor mínimo en la columna “xerror” de la tabla de complejidad cptable del árbol de decisión. Luego, con ese índice, se extrae el valor correspondiente de la columna “CP” en la tabla de complejidad.

plotcp(arbol)

Podamos el árbol como indica la función.

pruneTREE1 = prune(arbol, cp = arbol$cptable[which.min(arbol$cptable[,
"xerror"]), "CP"])
fancyRpartPlot(pruneTREE1, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")

pruneTREE2 = prune(arbol, cp = arbol$cptable[5, "CP"])
fancyRpartPlot(pruneTREE2, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")

tab2 = table(pred = predict(pruneTREE2, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
acierto2 = sum(diag(tab2))/ntest
acierto2
## [1] 0.8396275

Tras la poda nos sale literalmente es mismo resultado.

Creamos la matriz de confusión.

prediccion_1 <- predict(arbol, newdata = mhp_train, type = "class")
confusionMatrix(prediccion_1, mhp_train[["precio_bin"]])
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    0
##          1 3246  619
##          0  617 3250
##                                           
##                Accuracy : 0.8401          
##                  95% CI : (0.8318, 0.8482)
##     No Information Rate : 0.5004          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6803          
##                                           
##  Mcnemar's Test P-Value : 0.9773          
##                                           
##             Sensitivity : 0.8403          
##             Specificity : 0.8400          
##          Pos Pred Value : 0.8398          
##          Neg Pred Value : 0.8404          
##              Prevalence : 0.4996          
##          Detection Rate : 0.4198          
##    Detection Prevalence : 0.4999          
##       Balanced Accuracy : 0.8401          
##                                           
##        'Positive' Class : 1               
## 

Los resultados en la matriz de confusión son los siguientes:

  • Verdaderos positivos (VP): 3246
  • Verdaderos negativos (VN): 3250
  • Falsos positivos (FP): 619
  • Falsos negativos (FN): 617

También obtenemos:

  • Accuracy (exactitud): 84.01% de predicciones correctas.
  • Kappa: 68.03%
  • Sensitivity (sensibilidad): 0.8403. También conocida como recall o tasa verdadera positiva, mide la proporción de casos positivos que se identificaron correctamente.
  • Specificity (especificidad): 0.8400. Mide la proporción de casos negativos que se identificaron correctamente. Un valor más alto es mejor.
  • Pos Pred Value (valor predictivo positivo): 0.8398. Mide la proporción de predicciones positivas que son realmente positivas. Un valor más alto es mejor.
  • Neg Pred Value (valor predictivo negativo): 0.8404. Mide la proporción de predicciones negativas que son realmente negativas. Un valor más alto es mejor.
  • Mcnemar’s Test P-Value: 0.9773. Esta prueba compara las diferencias entre los falsos positivos y los falsos negativos. Un valor alto de “p” (como en este caso) sugiere que no hay una diferencia significativa entre ellos.

En resumen, este modelo de árbol de decisión podado tiene una exactitud del 84.01% y un coeficiente kappa de 0.6803, lo que indica un buen rendimiento en la clasificación.

Sin embargo para este modelo se han tenido en cuenta bastantes variables, y aunque el rendimiento es bueno, se puede alcanzar un rendimiento casi igual utilizando sólo los metros cuadrados, como se puede ver a continuación

arbol <- rpart(precio_bin ~ m2, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")

tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1
## [1] 0.8344542

Esta es una forma simple de confirmar la relevancia que tiene la superficie de una casa sobre su precio final, en contraste con otras variables como el tener ascensor o garaje. Además, esto nos deja un modelo mucho más simple y explicable.

10.4 Random Forest

Algoritmo de aprendizaje supervisado que se basa en la construcción de múltiples árboles de decisión durante el entrenamiento y la combinación de sus resultados para hacer una predicción. El objetivo principal de este enfoque es mejorar la precisión y la robustez del modelo en comparación con un único árbol de decisión.

El funcionamiento de un algoritmo Random Forest sería el siguiente:

  • Seleccionar muestras de arranque: El algoritmo selecciona aleatoriamente un subconjunto de muestras del conjunto de datos de entrenamiento, con reemplazo. Esto significa que una muestra puede ser seleccionada más de una vez.

  • Crear un árbol de decisión: Se crea un árbol de decisión utilizando el subconjunto de muestras seleccionado. Durante la construcción del árbol, en cada división del nodo, se selecciona un subconjunto aleatorio de características como candidatas para la división. La mejor división se elige en función de la medida de impureza como la entropía o el índice de Gini. Esto introduce aleatoriedad adicional en el proceso, lo que ayuda a reducir la correlación entre los árboles.

  • Repetir: El proceso se repite varias veces (número de árboles especificado) para construir un conjunto de árboles de decisión.

  • Combinar las predicciones: Para hacer una predicción en un nuevo caso, el algoritmo Random Forest ejecuta todos los árboles individualmente y luego combina sus resultados. En el caso de la clasificación, la clase con más votos (predicciones) de todos los árboles se elige como la predicción final.

Pasamos a construir el modelo comenzando con 10 árboles aleatórios.

# Escalamos la columna 4 en train y test
set.seed(19042023)
mhp_train_scaled <- mhp_train
mhp_train_scaled[, 4] <- scale(mhp_train[, 4])

mhp_test_scaled <- mhp_test
mhp_test_scaled[, 4] <- scale(mhp_test[, 4])

# Nos quedamos con todas las variables menos barrio y distrito porque hay muchas categorias
# Cogemos solo 10 árboles en el parámetro ntree
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
                           y = mhp_train_scaled$precio_bin,
                           ntree = 10)
# Predicción de los resultados con el conjunto de testing
y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])
confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    0
##          1 1676  127
##          0  257 1806
##                                           
##                Accuracy : 0.9007          
##                  95% CI : (0.8908, 0.9099)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8013          
##                                           
##  Mcnemar's Test P-Value : 4.61e-11        
##                                           
##             Sensitivity : 0.8670          
##             Specificity : 0.9343          
##          Pos Pred Value : 0.9296          
##          Neg Pred Value : 0.8754          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4335          
##    Detection Prevalence : 0.4664          
##       Balanced Accuracy : 0.9007          
##                                           
##        'Positive' Class : 1               
## 

Ahora probamos con 100 árboles aleatorios para ver si cambia mucho el resultado.

set.seed(19042022)
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
                           y = mhp_train_scaled$precio_bin,
                           ntree = 100)
y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])
confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    0
##          1 1684  108
##          0  249 1825
##                                           
##                Accuracy : 0.9077          
##                  95% CI : (0.8981, 0.9166)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8153          
##                                           
##  Mcnemar's Test P-Value : 1.267e-13       
##                                           
##             Sensitivity : 0.8712          
##             Specificity : 0.9441          
##          Pos Pred Value : 0.9397          
##          Neg Pred Value : 0.8799          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4356          
##    Detection Prevalence : 0.4635          
##       Balanced Accuracy : 0.9077          
##                                           
##        'Positive' Class : 1               
## 

La precisión del modelo es del 0,9077 y el intervalo de confianza del 95% es (0,8908, 0,9099), lo que indica que podemos estar bastante seguros de que la precisión real del modelo está dentro de este rango. El valor de kappa de 0,8013 indica que hay una buena concordancia entre las predicciones del modelo y las observaciones reales.

Además, la sensibilidad del modelo es del 0,8670 y la especificidad del 0,9343. La tasa de detección es del 0,4335 y la prevalencia es del 0,5000. En general, el modelo tiene un desempeño bastante bueno en la clasificación de las viviendas por encima y por debajo del umbral.

10.5 SVM

Support Vector Machine, algoritmo de aprendizaje supervisado que se utiliza tanto para problemas de clasificación como de regresión. En el caso de la clasificación, el objetivo principal de SVM es encontrar el hiperplano óptimo que separa los datos en diferentes clases. Un hiperplano es un límite de decisión que divide el espacio de características en dos partes, de modo que cada parte contenga puntos de datos de una clase específica. En el caso de la regresión, SVM busca encontrar una función que tenga el menor error de ajuste posible.

El algoritmo SVM funciona de la siguiente manera:

  • Encuentra el hiperplano óptimo: El algoritmo busca el hiperplano que tenga el margen máximo entre las dos clases. El margen es la distancia entre el hiperplano y los puntos de datos más cercanos a él, llamados “vectores de soporte”. Estos vectores de soporte son los puntos de datos que definen el límite de decisión y son fundamentales para el proceso de SVM.

  • Solución de un problema de optimización: La búsqueda del hiperplano óptimo implica la solución de un problema de optimización convexa. Este problema se resuelve mediante técnicas de optimización como el método de Lagrange o el algoritmo de optimización secuencial mínima (SMO).

  • Transformación a través del kernel: En casos en los que los datos no son linealmente separables en el espacio de características original, SVM puede utilizar una función de kernel para transformar los datos en un espacio de características de mayor dimensión donde se puedan separar linealmente. Los kernels comunes incluyen el kernel lineal, el kernel polinómico, el kernel de base radial (RBF) y el kernel sigmoide. La elección del kernel adecuado es crucial para el rendimiento del modelo SVM.

  • Clasificación: Para hacer una predicción en un nuevo punto de datos, el algoritmo SVM evalúa la posición del punto con respecto al hiperplano óptimo. Si el punto está por encima del hiperplano, se asigna a una clase, y si está por debajo, se asigna a la otra clase.

Para resolver de SVM necesito conocer el kernel y el parámetro de regularización que es una cuota superior sobre las variables duales, en definitiva, un coste. Pasamos a construir un primer el modelo usando un kernel radial RBF.

# Estandarizamos las variables numéricas "m2" y "rooms"
mhp_train$m2_estandarizado <- scale(mhp_train$m2)
mhp_train$rooms_estandarizado <- scale(mhp_train$rooms)

mhp_test$m2_estandarizado <- scale(mhp_test$m2)
mhp_test$rooms_estandarizado <- scale(mhp_test$rooms)

# Entrenamos el modelo SVM  kernel con kernel radial

modelo_svm <- svm(precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito, data = mhp_train, kernel = "radial", cost = 10, scale = TRUE)

# Hacemos predicciones en el conjunto de prueba
svm_predicciones <- predict(modelo_svm, mhp_test[, c("rooms_estandarizado", "m2_estandarizado", "garage", "elevator", "tipo_casa", "precio_distrito")])


# Calculamos la matriz de confusión y la precisión del modelo
svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test$precio_bin)
print(svm_matriz_confusion)
##           Real
## Prediccion    1    0
##          1 1692  116
##          0  241 1817
precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
print(paste0("Precisión del modelo SVM: ", round(precision * 100, 2), "%"))
## [1] "Precisión del modelo SVM: 90.77%"

Esto significa que:

  • 1692 casas con precio alto fueron clasificadas correctamente como precio alto (verdaderos positivos).
  • 1817 casas con precio bajo fueron clasificadas correctamente como precio bajo (verdaderos negativos).
  • 241 casas con precio alto fueron clasificadas incorrectamente como precio bajo (falsos negativos).
  • 116 casas con precio bajo fueron clasificadas incorrectamente como precio alto (falsos positivos).

La precisión del modelo es del 90.77%.


ELECCION DEL KERNEL

La elección del kernel en un modelo SVM depende de la naturaleza y la distribución de los datos. No hay una regla única para seleccionar el kernel, pero aquí hay algunas pautas generales para ayudarte a tomar una decisión:

  • Kernel lineal: Si tus datos son linealmente separables o la cantidad de atributos es alta en comparación con la cantidad de muestras, entonces un kernel lineal suele funcionar bien. Además, es computacionalmente más eficiente que otros kernels.

  • Kernel polinomial: Si tus datos tienen una estructura que puede ser capturada por un polinomio de un grado específico, el kernel polinomial podría ser una buena opción. Sin embargo, a medida que aumenta el grado del polinomio, el riesgo de sobreajuste aumenta y la complejidad computacional puede crecer.

  • Kernel radial (RBF): El kernel radial es una opción popular y versátil en muchos casos, especialmente cuando no se sabe de antemano si los datos son linealmente separables. Puede manejar datos no lineales y complejos. Sin embargo, encontrar los parámetros adecuados (costo y gamma) puede requerir una búsqueda en el espacio de parámetros y, en consecuencia, un mayor tiempo de cómputo.

  • Kernel sigmoide: El kernel sigmoide es similar a una función de activación de una red neuronal y puede ser una opción si deseas un enfoque similar al de una red neuronal pero con un modelo SVM.

Vamos a probar a representar los datos para ver si de alguna manera puede tomarse una mejor decisión sobre el Kernel a usar.

mhp_train %>%
  ggplot(aes(x = m2, y = rooms, color = as.factor(precio_bin))) +
  geom_point() +
  facet_grid(elevator ~ garage, labeller = labeller(
    elevator = c('0' = 'Sin Elevador', '1' = 'Con Elevador'),
    garage = c('0' = 'Sin Garaje', '1' = 'Con Garaje')
  )) +
  labs(title = "Diagrama de dispersión de m2 vs rooms, según garage y elevator",
       x = "Metros cuadrados",
       y = "Número de habitaciones",
       color = "Precio bin") +
  scale_color_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.title.position = "plot")


Probamos varios KERNELS

Creamos una función entrenar_evaluar_svm() que entrenará y evaluará un modelo SVM utilizando el tipo de kernel especificado en el argumento kernel_type. Luego con lapply() se aplicará esta función a cada tipo de kernel en el vector kernels. La función presentará un matriz de confusión y la precisión para cada tipo de kernel.

# Definir una función para entrenar y evaluar un modelo SVM con un tipo de kernel específico
entrenar_evaluar_svm <- function(kernel_type) {
  # Entrenamos el modelo SVM con el kernel especificado
  modelo_svm <- svm(precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito, data = mhp_train, kernel = kernel_type, cost = 10, scale = TRUE)

  # Hacemos predicciones en el conjunto de prueba
  svm_predicciones <- predict(modelo_svm, mhp_test[, c("rooms_estandarizado", "m2_estandarizado", "garage", "elevator", "tipo_casa", "precio_distrito")])

  # Calculamos la matriz de confusión y la precisión del modelo
  svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test$precio_bin)
  precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
  
  # Devolvemos una lista con la precisión y la matriz de confusión
  return(list(precision = precision, matriz_confusion = svm_matriz_confusion))
}

# Probar diferentes tipos de kernel
kernels <- c("linear", "polynomial", "radial", "sigmoid")
resultados <- lapply(kernels, entrenar_evaluar_svm)

# Imprimir los resultados
for (i in 1:length(kernels)) {
  cat(paste("Precisión del modelo SVM con kernel", kernels[i], ":", round(resultados[[i]]$precision * 100, 2), "%\n"))
  print(resultados[[i]]$matriz_confusion)
}
## Precisión del modelo SVM con kernel linear : 90.46 %
##           Real
## Prediccion    1    0
##          1 1706  142
##          0  227 1791
## Precisión del modelo SVM con kernel polynomial : 90.43 %
##           Real
## Prediccion    1    0
##          1 1688  125
##          0  245 1808
## Precisión del modelo SVM con kernel radial : 90.77 %
##           Real
## Prediccion    1    0
##          1 1692  116
##          0  241 1817
## Precisión del modelo SVM con kernel sigmoid : 85 %
##           Real
## Prediccion    1    0
##          1 1633  280
##          0  300 1653

Se observa que el KERNEL RADIAL es el que mejor predice, ahora vamos a ajustar algunos parámetros del modelo utilizando la función tune para encontrar los mejores parámetros.

# Entrenar un modelo SVM con diferentes parámetros
modelo_svm_radial <- svm(precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito, data = mhp_train, kernel = "radial", cost = 10, scale = TRUE)

# Utilizar la función tune() para encontrar los mejores parámetros
tuned_parameters <- tune(svm, 
                         precio_bin ~ rooms_estandarizado + m2_estandarizado + garage + elevator + tipo_casa + precio_distrito,
                         data = mhp_train,
                         kernel = "radial",
                         ranges = list(cost = c(0.1, 1, 10, 100, 1000),
                                       gamma = c(0.1, 0.5, 1, 2, 5)))

# Entrenar un modelo SVM con los mejores parámetros encontrados
best_model <- tuned_parameters$best.model

# Hacer predicciones en el conjunto de prueba con el mejor modelo
best_predictions <- predict(best_model, mhp_test[, c("rooms_estandarizado", "m2_estandarizado", "garage", "elevator", "tipo_casa", "precio_distrito")])

# Calcular la matriz de confusión y la precisión del mejor modelo
best_matriz_confusion <- table(Prediccion = best_predictions, Real = mhp_test$precio_bin)
print(best_matriz_confusion)
##           Real
## Prediccion    1    0
##          1 1692  118
##          0  241 1815
best_precision <- sum(diag(best_matriz_confusion)) / sum(best_matriz_confusion)
print(paste0("Precisión del mejor modelo SVM: ", round(best_precision * 100, 2), "%"))
## [1] "Precisión del mejor modelo SVM: 90.71%"

La precisión del modelo no mejora del 90.77% anterior y el coste computacional en la busqueda del mejor modelo es enorme. Parece que los parámetros del modelo elegidos durante el proceso de afinación no han sido los correctos. La afinación e hiperparámetros implica buscar la mejor combinación de hiperparámetros dentro del rango de búsqueda que se proporciona. Si el rango de búsqueda no contiene la combinación óptima de hiperparámetros, el proceso de afinación podría terminar eligiendo una combinación subóptima que resulte en un rendimiento inferior.

11 Evaluación y comparación de modelos

Tras la realización de todos los modelos, juntamos los mejores resultados en una sóla tabla para mejor comparabilidad:

11.0.1 K-fold Cross Validation

11.0.2 N-fold Cross Validation

12 Elección punto de corte